function SigSegmData = PluVelScalogram(sigPluVel,ParamWN,chPlot)
% PluVelScalogram - Scalogram computation for rainfall and velocity time series
%
% SigSegmData = PluVelScalogram(sigPluVel,ParamWN,chPlot)
%
% Let sigPluVel be either a time series of rainfall and velocity data or
% a time series of rainfall data only. In particular,  
% - In the first case it is sigPluVel = [t Plu Vel], where Plu(k) is the 
%   pluviometric value (in mm) corresponding to t(k), time espressed in 
%   MATLAB serial number form, with t(k+1)-t(k)=1/24, i.e. 1h, for each 
%   k, and Vel(k) is the corresponding velocity (the sampling frequency is 
%   Fs=24 cpd). Another requirement is that t(1) and t(end) are integer 
%   numbers, i.e. are dates related to 0h.
% - In the second case it is sigPluVel = [t Plu], with the same requirements
%   about t and Plu.
% The elements of the input Plu can be differential or cumulative values 
% (at some days or fraction of day), according to the user's preferences. 
%
% This function computes the scalograms by means of Continuous Wavelet 
% Transform (CWT). The parameters for CWT computation are taken from
% ParamWN, managed by means of DefParam (if ParamWN is undefined or empty, 
% DefParam is automatically called).
% For computation efficiency purposes, if 'amor' is the selected wavelet,
% the CWT is computed by means of preliminary CWT filter bank generation.
% If SigPluVel is a 3-columns matrix, each saved image (see below) shows 
% two scalograms: the upper one is related to rainfall data, and the lower 
% one is related to velocity data. If Nb > 0 (see below), the velocity 
% scalogram is zero-padded in order to align it to the rainfall scalogram. 
% If SigPluVel is a 2-columns matrix, each output image shows a scalogram 
% only, related to rainfall data.
%
% Each scalogram is computed for the reference date t(k) on the basis of 
% these input data taken from ParamWN:
%
%   - Na        : the starting date of the scalograms related to t(k) is
%                 t(k)-Na, where Na is expressed in days, under the condition 
%                 t(k)-Na >= t(1) (if this condition is not fulfilled, the 
%                 corresponding scalogram is not computed). 
%                 If Na is undefined or empty, the default value Na = 15 is
%                 used.
%
%   - Nb        : the ending date of the pluviometric scalogram related to 
%                 t(k) is t(k)+Nb, where Nb is expressed in days, under the 
%                 condition that t(k)+Nb <= t(end) (if such a condition is 
%                 not fulfilled, the corresponding scalogram is not computed).
%                 The ending date of the velocity scalogram always is t(k). 
%                 If Nb is undefined or empty, the default choice Nb = 0
%                 is used.
%
%   - NPerDay   : number of scalograms per day. The allowed values are
%                 1 (one scalogram per day at 0h), 2 (two scalograms at 
%                 0h and 12 h), 3 (0h, 8h, 16h), 4 (0h, 6h, 12h, 18h), 6
%                 (0h, 4h, 8h, 12h, 16h, 18h, 20h) 12 (even hours), 24.    
%                 If nPerDay is undefined, empty or not an allowed value, 
%                 nPerDay = 1 is used.  
%
%   - DateIn    : initial date for which the scalograms are computed (this
%                 date must be expressed as serial number). The start time 
%                 of the first computed scalogram is max(d1-Na,t(1)). 
%                 If d1 is undefined of empty, d1 = t(1).
%
%   - DateFin   : final date for which the scalograms are computed (this
%                 date must be expressed as serial number). The end time 
%                 of the last computed scalogram is min(d2+Nb,t(end)). 
%                 If d2 is undefined of empty, d2 = t(end).
%
%   - SizeIm    : size of generated image files. 
%
%   - ComPart:    common part of the output filename (folder name, constant 
%                 part of the filename). Each filename is completed with the 
%                 date string. For example, if comPart is the string 
%                       'pluvioData\LongBeach15days', 
%                 the filename of the scalogram of data related to the Na 
%                 days until the date 15 march 2020, 12:00:00 is the string
%                    'pluvioData\LongBeach15days_15-Mar-2020-120000.jpg'.
%                 If comPart is undefined or empty, comPart = '' is used.
%
% If ParamWN is a char variable, it is supposed to be the name of a file
% having a field named 'ParamWN'. If the file loading or the extraction
% of the ParamWN value is not successful, it is ParamWN = [].
% If ParamWN is undefined or empty, a session of DefParam runs to generate
% ParamWN.
%
% The remaing input argument is:
%
%   - chPlot:     vector of dates for which, besides the file generation, 
%                 the scalogram are also shown in a two subplots figure,
%                 where the upper scalogram is the rainfall one and the 
%                 lower scalogram is the velocity one (if rainfall data
%                 only are used, a scalogram only is shown). 
%                 The dates in chPlot must be coherent with the choice on
%                 NPerDay.
%                 If chPlot is undefined or empty, no figures are shown.
%
% The duration of each rainfall signal segment is Na+Nb days and the 
% corresponding length is (Na+Nb)*Fs+1. The duration of each valocity 
% signal segment is Na days, with Na*Fs+1 length.
%
% The output SigSegmData is a 4-column cell array such that SigSegmData{k,1} 
% is is the reference time of the k-th row, SigSegmData{k,2} is the initial
% time of the rainfall time series, SigSegmData{k,3} is the corresponding
% end time, SigSegmData{k,4} is the generated filename. In this way, some
% information about the signal segment which contributes to the k-th figure 
% can be stored.
%
% See also PluVelInspect, trendClass, dataHomAug, dataTVT, TRLearnPluVel,
%   DefParam.
%
% SigSegmData = PluVelScalogram(sigPluVel,ParamWN,chPlot)

% G. Teza, 2020

if nargin < 3
    chPlot = [];
end

if nargin < 2
    ParamWN = [];
end
if ischar(ParamWN)
    try 
        sv = load(ParamWN);
        ParamWN = sv.ParamWN;
    catch
        ParamWN = [];
    end
end
if isempty(ParamWN)
    ParamWN = DefParam;
end

comPart = ParamWN.ComPart;
[pathFold,~,~] = fileparts(comPart);
if ~isempty(pathFold)
    if exist(pathFold,'dir') ~= 7
        mkdir(pathFold)
    end
end

Wavelet = ParamWN.Wavelet;

sizeIm = ParamWN.SizeIm;

d1 = ParamWN.DateIn;
d2 = ParamWN.DateFin;

Na = ParamWN.Na;
Nb = ParamWN.Nb;

nPerDay = ParamWN.NPerDay;
Fs = ParamWN.Fs;

t = sigPluVel(:,1);
dt = median(diff(t));
if abs(dt-1/Fs) > 10^6*eps
    error('Invalid input data: the sampling rate must be 1/%d cpd',Fs);
end

sigP = sigPluVel(:,2);
lSigP = round((Na+Nb)*Fs)+1;    % length of each pluvio segment
if size(sigPluVel,2) == 3
    sigV = sigPluVel(:,3);
    lSigV = round(Na*Fs)+1;     % length of each velocity segment
else
    sigV = [];
    lSigV = 0;
end
if isempty(d1) || (d1 < t(1))
    n1 = 1;
else
    d1 = floor(d1);
    [~,n1] = min(abs(t-d1));
end
if isempty(d2) || (d2 > t(end))
    n2 = numel(t);
else
    d2 = ceil(d2);
    [~,n2] = min(abs(t-d2));
end

nStep = round(24/nPerDay);

% Filter bank (applicable only if default 'amor' is used)
if strcmp(Wavelet,'amor')
    % CWT filter bank for pluvio data
    cwtfbP = cwtfilterbank(...
        'SignalLength',lSigP,...
        'Wavelet','amor',...
        'VoicesPerOctave',24,...
        'samplingFrequency',Fs);
    % CWT filter bank for velocity data (if applicable)
    if ~isempty(sigV)
        cwtfbV = cwtfilterbank(...
            'SignalLength',lSigV,...
            'Wavelet','amor',...
            'VoicesPerOctave',24,...
            'samplingFrequency',Fs);
    end
end

nt = n2-n1+1;

hwb = waitbar(0,'Scalogram generation in progress...');
cmap = jet(128);
nplot = 1;      % for the scalogram plot, if applicable
SigSegmData = cell(numel(n1:nStep:n2),4);
ncont = 1;      % for the correct indexing of SigSegmData
for k = n1:nStep:n2
    trefk = t(k);
    t1k  = trefk-Na;
    t2Pk = trefk+Nb;
    if (t1k >= t(1))&&(t2Pk <= t(end))
        [~,n1k]  = min(abs(t-t1k));
        [~,n2Pk] = min(abs(t-t2Pk));
        tPk = t(n1k:n2Pk);
        sigPk = sigP(n1k:n2Pk);
        if ~isempty(sigV)
            t2Vk = trefk;
            [~,n2Vk] = min(abs(t-t2Vk));
            tVk = t(n1k:n2Vk);
            sigVk = sigV(n1k:n2Vk);
        else
            sigVk = [];
        end
        if nonanTS(sigPk) && nonanTS(sigVk) && numel(sigPk) == lSigP && numel(sigVk) == lSigV
            if strcmp(Wavelet,'amor') % Computation by using the precomputed filter banks:
                [cfsPk,fPk] = cwtfbP.wt(sigPk);
                if ~isempty(sigV)
                    [cfsVk,fVk] = cwtfbV.wt(sigVk);
                end
            else
                Omega0 = 2*pi*F0;
                [cfsPk,tPk] = contwt(sigPk,1/Fs,[],0.13,[],[],'morlet',[Omega0 Sigmat]);
                fPk = 1./tPk;
                if ~isempty(sigV)
                    [cfsVk,tVk] = contwt(sigVk,1/Fs,[],0.13,[],[],'morlet',[Omega0 Sigmat]);
                    fVk = 1./tVk;
                end
            end
            % Plot (if required):
            if ~isempty(chPlot) && (nplot <= numel(chPlot))
                if abs(t2k-chPlot(nplot)) < 10^3*eps
                    %----------------------------------------------
                    plotScalograms(tPk,fPk,cfsPk,tVk,fVk,cfsVk);
                    %----------------------------------------------
                    nplot = nplot+1;
                end
            end
            % Image generation ad saving:
            filenak = [comPart '_' datestr(trefk) '.jpg'];
            filenak = strrep(filenak,' ','-');  % to change " " to "-" in datestr  
            filenak = strrep(filenak,':','');   % to remove ":"  in datestr
            imPk = ind2rgb(im2uint8(rescale(abs(cfsPk))),cmap);
            if ~isempty(sigV)
                imVk = ind2rgb(im2uint8(rescale(abs(cfsVk))),cmap);
                if Nb > 0   % zero padding of velocity scalogram, if necessary
                    MZeroPad = zeros(size(imVk,1),size(imPk,2)-size(imVk,2),3);
                    imVk = cat(2,imVk,MZeroPad);
                end
                imk = cat(1,imPk,imVk);
            else
                imk = imPk;
            end
            imk = imresize(imk,sizeIm);
            imwrite(imk,filenak);
            waitbar((k-n1)/nt,hwb);
            SigSegmData(ncont,:) = {trefk,t1k,t2Pk,filenak};
            ncont = ncont+1;
        end
    end
end
close(hwb);

if ncont < size(SigSegmData,1)      % to remove possible useless rows
    SigSegmData(ncont:end,:) = [];
end


%% Ancillary functions:

function Ina = nonanTS(TS)
% Ina is 1 if TS does not contain any NaN or Inf value, 0 elsewhere. 
I1 = isnan(TS);
I2 = isinf(TS);
if (sum(I1) > 0) || (sum(I2) > 0)
    Ina = 0;
else
    Ina = 1;
end
Ina = logical(Ina);